Temperature dependence of electrokinetic flux in Si nanochannel 
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Significant temperature effects on the electrokinetic transport in a nanochannel with a slab 
geometry are demonstrated using a molecular dynamics (MD) model. A system consisting of Na + 
and CP ions dissolved in water and confined between fixed crystalline silicon walls with negatively 
charged inner surfaces in an external electric field was investigated. Lennard- Jones (LJ) force fields 
and Coulomb electrostatic interactions with Simple Point Charge Extended (SPC/E) model were used 
to represent the interactions between ions, water molecules, and channel wall atoms. Dependence of 
the flow of water and ions on the temperature was examined. The magnitude of the water flux and 
even its direction are shown to be significantly affected by temperature. In particular, the previously 
reported flow reversal phenomenon does not occur at higher temperature. Temperature dependence 
of the flux was attributed to the charge redistribution and to the changes in viscosity of water. 
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I. INTRODUCTION 



Nanoscale numerical models of electro-osmosis [THS] 
provide insight into the important transport mechanism, 
help improve current technological devices, and guide the 
design of new technology based on the principles of elec- 
trokinetic transport. They also improve understanding of 
transport processes in biology [7] and their temperature 
dependence [5]. 

The present study examines the fluid flow and ion 
species transport under an electric field in nanochannels 
typically found in heterogeneous porous media. Although 
the continuum conservation equations can be applied to 
microscale channels, the main driving force of electroki- 
netic transport occurs in an electric double layer (EDL) 
at the solid liquid interface with dimensions that can be 
comparable to intermolecular distances. Therefore, molec- 
ular dynamics (MD) simulations were applied to analyze 
the interaction between ions, water molecules, and wall 
atoms in the EDL region. Time averaged velocity and 
concentration profiles of water molecules and ionic species 
at different temperatures were obtained from these MD 
simulations. Water viscosity profile across the channel 
was estimated at different temperatures and its impact on 
the transport of water and ionic species will be discussed. 




|Video i"j (Color online) Simulation box without water 
molecules. Si wall atoms are gray, Na + ions red, and Cl~ 
ions green. Water molecules are not shown. Accompany- 
ing animation shows movement of ions when electric field is 
applied in the positive x direction. 



II. MD POTENTIALS 

Models based on classical Lennard- Jones (LJ) force 
fields and Coulomb electrostatic interactions with Simple 
Point Charge Extended (SPC/E) [9] model were used to 
represent the interactions between ions, water molecules, 
and channel wall atoms. LJ contribution of atoms i and 
j to the total potential energy is 



V LJ {r) 
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The parameters (the depth of the potential well) 
and &ij {\f2oij is the minimum energy distance) depend 
on atomic species of ith and jth atoms (Table [T]). The 
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TABLE I. Parameters of Lennard-Jones potentials (from the 
GROMACS [11] force field), ov, in A, ey in cal/mol. 



ij 


oo 


OSi 


ONa 


OC1 


SiSi 


SiNa 


SiCl 


NaNa 


NaCl 


C1C1 


CTij 
Eij 


3.17 
155 


3.27 
301 


2.86 
47.9 


3.75 
129 


3.39 
584 


2.95 
92.9 


3.88 
249 


2.58 
14.8 


3.38 
39.6 


4.45 
106 



Coulomb electrostatic potential energy contribution is 

1 QtQj 



V c {r) 



47T£n 



(2) 



The £q is vacuum permittivity, q^ and qj are charges of ith 
and jth atoms, and r is the distance between ith and jth 
atoms. The Particle-Particle Particle-Mesh (PPPM) [10] 
method was used for long range electrostatics. 



III. SIMULATION SETUP 

In order to study the temperature effects, the authors 
first reproduced the system studied previously [2], a sum- 
mary of which is given next. The dimensions of the 
solution region were 4.66x4.22x3.49 nm. Channel walls, 
perpendicular to the z axis, were formed by four [111] 
oriented layers of Si atoms in a diamond crystal structure, 
each wall being 0.39 nm thick. Periodic boundary con- 
ditions were applied in the x and y directions. The size 
of the simulation cell in the z direction was extended to 
three times outermost-to-outermost wall layer distance 
(3x4.37 nm) to mitigate electrostatic interactions of pe- 
riodic images in the z direction. The electric field of 
0.55 V/nm was applied in the positive x direction. 

The crystalline channel walls consisted of 1232 fixed 
Si atoms. A charge of —70 e was distributed uniformly 
on atoms of the innermost surface layers of both Si walls 
(total 308 Si atoms were charged negatively and 924 were 
uncharged). Therefore, the innermost wall surface atoms 
had —0.227273 e/atom charge. That corresponds to a 
surface charge density of —0.285 C/m 2 , which is close 
to typical value for a fully ionized surface 0.2 C/m 2 [12] 
or the charge density of —0.1 C/m 2 measured at silica 
surface [T3]. Then 2290 water molecules were inserted 
avoiding close contacts, and randomly selected 146 water 
molecules were replaced by 108 Na + and 38 CD ions. 
Thus the whole system was electrically neutral and re- 
sulting NaCl concentration in the channel center at the 
temperature of 300 K was sal. 2 M. 

First, the energy minimization of the system was per- 
formed using the conjugate gradient method. Then, the 
system was equilibrated by 2 ns of MD simulation without 
an electric field. A timestep of 2xl0 -15 s was used for the 
leapfrog [14] integration of Newton's equations of motion. 
The resulting configuration, excluding water molecules, is 
shown in Video [T] Finally, a 22 ns MD run was performed 
with external electric field. The SHAKE [15 algorithm 
was used to constrain bonds of water molecules. The 
solution temperature was controlled by the Nose-Hoover 



IV. THERMOSTAT 

The thermostat in MD simulations adjusts the veloc- 
ity to maintain the desired temperature and its proper 
fluctuations. To check how much the velocity compo- 
nent in the field direction is affected by thermostat, the 
velocity profiles from simulation with selective thermo- 
stat (adjusting y and z velocity components only) were 
compared to results with full thermostat (adjusting x, y, 
and z components of velocity) . The number of thermal 
degrees of freedom was adjusted accordingly. Figure [I] 
shows that selective thermostat produces slightly more 
negative velocity. 

Velocity profiles with profile-unbiased thermostat 
(PUT) [TT] that "preserves" the velocity profile (i.e. its 
x component) across the channel were also generated. 
Averages did not differ significantly from those obtained 
by selective thermostat (adjusting only y and z velocity 
components), but PUT error bars were larger, therefore 
the selective thermostat (adjusting y and z components 
of velocity only) was used. 

The authors then verified that velocity profiles with full 
thermostat obtained from GROMACS [TT] were statistically 
identical to those from lammps [18] when thermostating 
all velocity components, even though GROMACS used the 
PME (Particle-Mesh-Ewald) method Q]5] with a slab cor- 
rection [3U] in the z direction for long range electrostatics 
and SETTLE [21] algorithm to constrain bonds of water 
molecules. 
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FIG. 1. Comparison of water velocity profiles using full 
thermostat (adjusting x, y, and z velocity components) and 
selective thermostat (adjusting y and z components only). 
Error bars are only shown on one side of the values to avoid 
overlap. The values shown (on all figures) are averages over 
~0.026 nm wide bins parallel to the xy plane that are then 
symmetrically averaged about the channel center. 

There is a noticeable difference between the velocity 
profile in Fig. [I] and the one of [2] ■ The velocities from 
the present work are more positive. That is consistent 
with the positive peak in the driving force 0.65 nm from 
the channel wall in Fig. [3j whereas the driving force of [5] 
remain negative in that region. Assuming that all the 
parameters were set correctly, the difference is hard to 
track since the gromacs [TT] package used by [2] does 
not offer a selective thermostat. 
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V. CONCENTRATION PROFILES 



VI. DRIVING FORCE 



The calculated atomic concentration profiles shown in 
Fig. [2] which agree well with reported results [5] , manifest 
formation of alternatively charged layers of atoms paral- 
lel to the channel walls. The negatively charged Si wall 
attracts both positively charged Na + ions and slightly pos- 
itively charged H atoms (0.4238 e/atom for SPC/E water) 
from H2O molecules, forming the first near-wall concen- 
tration peak. The adjacent layer is formed of slightly 
negatively charged O atoms (—0.8476 e/atom for SPC/E 
water) from water molecules. Five noticeable layers with 
alternating charge signs are formed — four ionic layers can 
be seen in Fig. [3] and the fifth is a negatively charged 
layer of O atoms 0.25 nm from the channel wall (Fig. [2]). 
Layering of particles near a flat surface is characteristic for 
polar liquids [22] , L J fluids [23] , and charged surfaces [24] . 
Further towards the channel center, the concentrations of 
ions are more balanced (sal. 2 M). 
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FIG. 2. Concentrations of Na + and Cl~ ions, and H and O 
atoms from H2O molecules symmetrically averaged across half 
of the channel. The inside plot is a magnified version of a 
larger plot. Concentration of Na + ions in the channel center 
increases as the temperature increases. 



Since the net charge of a standalone water molecule is 
zero, its center of mass will not move in the presence of 
external electric field. The driving force for the electro- 
osmotic flow comes from the electric field causing the 
movement of a fluid region with non-zero net charge. 
Regions of fluid with positive net charge will drive the 
flow in the direction of an external electric field, while 
the regions with negative net charge will drive the flow in 
the direction opposite to an external electric field. The 
driving force is defined as 



F d (z) = e [c Na + (z) - cci- (z)} E e 



(3) 



where e is the elementary charge, c Na + (z) and c cl - (z) are 
ionic number densities across the channel, and E ext is an 
external electric field. Figure [3] shows a dependence of 
the driving force on temperature. 

When the temperature is increased, the driving force 
in the region further than 0.42 nm from the channel 
wall becomes more positive (because of increased Na + 
concentration in that region), and water starts to flow 
in the positive x direction. The driving force from Cl _ 
ions, on the contrary, will increase only in the near-wall 
region, as some of the Cl~ ions will redistribute closer 
than 0.5 nm within the channel wall. In the following 
section it will be argued that this charge redistribution 
at increased temperature drives the flow in the positive x 
direction. 
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FIG 3. Driving force for the fluid flow given by Eq. (pjjl. 
Values are symmetrically averaged across half of the chan- 
nel. Inside plot magnifies values in the central region of the 
channel — where driving force becomes more positive at higher 
temperature. The temperature changes in concentration pro- 
files, resulting in changes of driving force, are the main factor 
contributing to the temperature dependence of the flow. 
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FIG. 4. Velocity profiles from MD at T=280 K (top), 300 K, and 320 K (bottom). Negative flow means movement in the 
direction opposite to the applied electric field. Velocity of H2O is a velocity of oxygen atoms from H2O molecules. Note that the 
temperature increase changes the water flow direction along the channel. Error bars (for Na + and H2O they are represented by 
rectangles) were obtained by analysis of ten simulations with different seeds of initial Gaussian random velocity distribution. 
Large error bars near the wall are due to low concentrations (Fig. [21) . Details of water velocity profiles are shown in Fig. [5] 



VII. TEMPERATURE DEPENDENCE OF 
VELOCITY PROFILES 

A significant dependence of the flow on temperature 
was observed. Figure |4] shows the dependence of the ve- 
locity profile of ions and water on the temperature. It 
demonstrates that the water flow and even its direction 
are directly affected by the temperature. At lower temper- 
atures, flow reversal of water was observed, in agreement 
with previous results [2J. However, at higher temper- 
atures, the water flows in the direction expected by a 
standard EDL model. 

Note that even though CI - ions in water move faster 
than Na + ions (Fig.[4|, that does not affect water velocity 
profile. Water movement is not governed by the ionic 
velocities, but by the profile of the driving force. Velocity 
of ions relative to water is given by their respective mobil- 
ities fii according to the definition of mobility Vi = /i^E. 
The LJ potential properly reflects the experimental mobil- 
ity of CI - ions being higher than the mobility of Na + ions. 
Mobility of ions will affect the rate of ion accumulation on 
electrodes in an experiment, but not the velocity profile 
of water in the present model. 

The temperature dependence of the electro-osmotic 



flow was attributed to (a) charge redistribution and to 
(b) changes of water viscosity with temperature. 



As mentioned in Sec. |VI[ the positive ions will redis- 
tribute towards the channel center at higher temperatures 
(Fig. [2}[3|. At T=320 K, Na+ ions dominate CT ions 
in the region 0.61 to 0.76 nm from the wall. The in- 
creased Na + concentration in the channel center comes 
from lowering two near-wall concentration peaks of Na + 
ions (at 0.14 and 0.38 nm from the wall in Fig. [2]), that 
redistribute further than 0.42 nm from the channel wall at 
higher temperature. This increased Na + concentration in 
the channel center (where the water viscosity is lower than 
at the surface 



see Sec. 



IX|) stimulates the movement of 

ions 



water in the positive x direction. Even though Na 
dominate only in regions less than 0.52 nm and 0.61 to 
0.76 nm from the wall, hydrogen bonding and collisions 
will also drag adjacent layers of water in the channel cen- 
ter, outperforming the competing mechanism of water 
dragged in the negative x direction by CI - ions that have 
higher concentration than the Na + in the channel center. 

On the contrary, at lower temperatures CI - ions are 
more dominant over Na + ions in the channel center, caus- 
ing a flow in the negative x direction. 
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VIII. RELATION OF VELOCITY PROFILES 
AND CHARGE DISTRIBUTION 

Velocity profile is related to charge density profile by 
the Stokes equation [25 



T=280K 



_d_ 

dz 



du x (z) 
dz 



= -F d (z). 



(4) 



The magnitude of the driving force, Fd(z), given by 
Eq. ([3]), as well as the velocity profile, u x (z), can be 
calculated as averages from MD simulations. The f](z), 
viscosity profile of water across the channel, can then be 
estimated from MD as follows. 



IX. ESTIMATION OF WATER VISCOSITY 

The water viscosity profile across the channel was es- 
timated following the method of [T] . The same method 
(except smoothing of velocity profile) was used by [21)] . 
To simplify equations, the coordinate system with z=0 
at the channel center will be used. Integration of Stokes 
Eq. Q leads to viscosity estimate 







- / F d (z)dz 


_ Jo 




du x {z) 




dz 


z=z 



(5) 



The fit diverges near the points where the derivative of 
velocity is zero. This is appropriate, since the viscosity 
can not be estimated in the region with zero shear strain — 
which is the term in denominator. The numerator term 
represents the shear stress [26] . 

To obtain smooth derivative of u x {z), the velocity pro- 
file was approximated by the sum of harmonic components 



Uxfit (z) = ^2 



j a n cos 

n=0 



(6) 



The h is the distance of the furthest oxygen atom from 
the channel center. In contrast to the original approxi- 
mation of Q], the exponential term is excluded. Also, to 
exploit the symmetry of the system, the origin of cosine 
components is set to the channel center. 

The velocity, its approximation, the viscosity profile 
estimated by Eq. ([5]), and the experimental viscosity of 
water [27 at the simulated temperatures are shown in 
Fig. [5j Note that the estimated viscosity (the line made 
of + symbols on the upper plots in Fig. [5]) at its minimum 
reproduces the experimental viscosity (solid horizontal 
line). 

At lower temperatures, the viscosity of near-wall water 
layers increases drastically. That means the near-wall 
Na + ions can not drive much water flow. Conversely, 
at higher temperatures, the near-wall water will become 
more mobile, exhibiting a partial slip at T=320 K. 
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FIG. 5. Velocity profile fit to the sum of harmonics (lower 
plots) and viscosity estimate from Eq. |5| (upper plots) at 
T=280 K, 300 K and 320 K. The inside plots are magnified 
versions of larger plots. 
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X. CONCLUSIONS 

Many technological processes and phenomena in live 
organisms exhibit a strong dependence on temperature. 
It is important to understand how the temperature affects 
technological and natural processes. This work indicates 
that a variation of temperature can have deterministic 
effects on an electrokinetic flow. 

This paper revisited the phenomenon of flow reversal 
during electrokinetic flow in a slit nanochannel. It was 
found, using a molecular dynamics analysis, that both the 
magnitude and direction of the electro-osmotic flow are 
significantly affected by temperature. Even though the 
results reported in [2] could be replicated, it was found 
that the flow reversal of water molecules no longer occurs 
if the temperature is increased slightly above the 300 K 
used in [2] . The mechanisms that lead to such a significant 
temperature dependence of the nanochannel flow were 
then investigated. In particular, it was shown that the 



temperature dependence of the flux can be explained by 
the charge redistribution and the decrease of near-wall 
water viscosity at higher temperatures. 
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